Weathering

Broadly, weathering is the process by which minerals are broken down. Weathering is an important part of many biogeochemical cycles and is responsible for many of the common constituents found in streamflow.

Chemical vs physical

There are two kinds of weathering and water is top level control on both.

Physical weathering is the mechanical breakdown of minerals via physical force. A simple example of physical weathering is runoff water picking up and carrying soil into a nearby river. Physical weathering can be carried out by any number of vectors: wind, water, freezing and thawing, tectonic shifts, etc. Watershed scientists sometimes use total suspended soils (TSS) or something similar as an imperfect proxy for this. The USGS is working on the next generation of suspended soil detection using acoustic methods. If you want to learn more you can read about it here.

Physical weathering is often the first steps in complete weathering of a mineral. With many of the later steps occurring chemically.

Chemical weathering is the reaction and dissolution of minerals. Water is a key player here because it acts as an acid, a base, and a medium for dissolution. Water can act as an acid, splitting chemical bonds or as a redox agent. Water plays a particularly key role in chemical weathering, the convection of dissolved materials. As you learned in the thermodynamics unit, once a reaction reaches equilibrium it has no further net change. Water constantly flowing through rock and soil carries away dissolved consitutents, keeping weathering reactions from ever reaching equilibrium. There are myriad ways minerals are weathered chemically and each can be sensitive to other environmental factors (more on that in a bit). In order to get at chemical weathering rates, watershed scientists will look at particular rock-derived solutes.

Rock derived solutes

Many common dissolved constituents of stream water are weathered from rock. Some of these, like calcium, sodium, potassium, manganese, and phosphorus, are nutrients needed to sustain life. Others, like aluminum, are toxic to life. All are found in their dissolved form due, primarily, to chemical weathering.

Most Abundant Element Approximate % by Weight
O 46.60
Si 27.70
Al 8.10
Ca 5.00
Na 3.70
K 2.60
Mg 1.50
Ti 0.44
P 0.10

Table of elemental composition of the Earth’s crust. Source: Wikipedia contributors. ‘Earth’s crust.’ Wikipedia, The Free Encyclopedia, 20 Feb. 2022. Web. 3 Mar. 2022.

Flux II: Return of the Flux

In our last R assignment (Gas Flux) we calculated the flux of methane off of the Yahara River. In this assignment we will be looking at a different kind of flux: solute flux from a watershed (often called ‘export’). Fluxes from terrestrial watersheds are snapshots a key linkage in global biogeochemical cycles, connecting land, freshwaters, and the ocean.

We can calculate solute flux much more easily than gas flux. This is because instead of grasping weakly at an estimate for transfer velocity (k) we just get to use discharge, something we are generally good at measuring. The equation for an instantaneous solute flux (like a snapshot in time) is:

\(Flux = (Q*[X])/A\)

Where ‘Q’ is discharge, ‘[X]’ is the concentration of a solute, and ‘A’ is the watershed area. Fluxes are usually reported in a mass per time per area.

Scientists and managers are often interested in the total flux over a period of time.The equation for solute flux over a given period is:

\(Flux = \sum_{t=1}^{t_x} (Q*[X])/A\)

Where ‘t’ is time. It’s important to consider the time step of an instantaneous flux before scaling it up to a larger time step.

Does solute flux == weathering?

When thinking about solute fluxes as a proxy for weathering rates we need to be careful. While many solutes are originally derived from rock, they may come from other sources on short enough timescales.

  • Anthropogenic: Humans can add solutes in many forms. For example, phosphorus is often used as a fertilizer and in de-icers for roadways.

  • Biologic: Solutes concentrations may also be affected by biological demand. Manganese is an important nutrient, meaning plants are actively trying to take it up.

  • Atmospheric: Dust and rain can deposit materials weathered elsewhere into the watershed. For example, sites near the coast may experience heavy sodium deposition from sea spray.

  • Complexation: As you learned in the Acid Rain assignment, some solutes (like aluminum) have solubilities that depend on the surrounding redox state or pH.

An ideal solute for estimating weathering rates will either not experience any of these factors or do so in a predictable way.

Importance over human timescales

Weathering, and the solute fluxes derived from it, are deeply important on the human (less than ~100 year) timescale. At it most basic level, weathering rates determine the baseline concentrations of many key solutes for human and aquatic health. They also represent a key indicator of soil health, showing land managers which elements are being used efficiently by plants and which are going to waste. Weathering rates also control the availability of phosphorus in natural systems, acting as an extended control on plant growth.

Importance over geologic timescales

Over geologic time (millions to billions of years) weathering plays a critical role in many biogeochemical processes. For example, calcium weathered from rock eventually reacts with dissolved bicarbonate to form calcium carbonate, the mineral sand and shells are made from. Remember from the acidity and alkalinity unit that bicarbonate is largely formed from dissolved, atmospheric carbon dioxide. This cycle of greater weathering leading to more buried carbon acts a stabilizing feedback loop for the earth climate system.

Conceptual diagram of stablizing weathering-climate feedback loop. Frings, Patrick J. ; Palaeoweathering: How Do Weathering Rates Vary with Climate?. Elements 2019;; 15 (4): 259–265. doi: https://doi.org/10.2138/gselements.15.4.259

Some researchers are looking into how we can use this effect over more relevant timescales to humanity. You can read more about it in the Taylor et al 2017 publication in the ‘pubs’ folder of this project.

Assignment

In this assignment, we are going to look at data from the Bear Brook Watershed in Maine (BBWM). The BBWM experiment has been running since the mid 1980s and is broadly focused on investigating how increased nitrogen, sulfur, and acidity affect surface waters and their surrounding watersheds.

In a previous assignment, we looked at how acid rain was changing over time and how it can control aluminum toxicity. The BBWM experiment is essentially testing what would happen if acid rain had continued unabated. BBWM is split into two watersheds, East Bear (EB) and West Bear (WB). EB is an untreated reference watershed. WB had ammonium sulfate dumped onto it on a bimonthly basis from 1989-2016. Ammonium sulfate is essentially a mixture of sulfuric acid (H2SO4) and ammonium (NH4+). This means there are a few things going on in WB:

You can see the watersheds here:

Q1

Part A

Look at the equation for solute flux calculation above. Why do we report fluxes adjusted by area?

Part B

Given the background information above and the controls on weathering covered in lecture, write a brief (1-2 paragraphs) hypothesis about how and why you think weathering rates will differ between EB and WB.

Make sure your answer covers both what you expect to differ and why you think that will be the case.

Part C

Consider the three affects listed in the background section of this assignment (increased nitrogen, sulfur, and acid inputs). For each, give an example of where you might find this condition in the real world. Acid rain from industrial NOx and SOx has been thoroughly covered already, so pick other examples.

Q2

Data loading and manipulation.

Part A

Read in the discharge (Q) and stream chemistry data for both watersheds (see the ‘data’ folder). Combine them into a single data frame and remove the quality assurance flags from the data. Use the function glimpse() to show your final, combined data frame.

(Hint: This is the same Macrosheds data format as the Acid Rain assignment. Go look at that if you are having trouble.)

data <- read_feather('data/EB_chem_stream.feather') %>%
  rbind(read_feather('data/WB_chem_stream.feather')) %>%
  rbind(read_feather('data/WB_q.feather')) %>%
  rbind(read_feather('data/EB_q.feather')) %>%
  select(datetime, site_code, var, val)

glimpse(data)
## Rows: 227,960
## Columns: 4
## $ datetime  <dttm> 1988-03-16, 1988-03-17, 1988-03-18, 1988-03-19, 1988-03-20,~
## $ site_code <chr> "EB", "EB", "EB", "EB", "EB", "EB", "EB", "EB", "EB", "EB", ~
## $ var       <chr> "GN_Al", "GN_Al", "GN_Al", "GN_Al", "GN_Al", "GN_Al", "GN_Al~
## $ val       <dbl> 0.076, 0.082, 0.088, 0.094, 0.101, 0.107, 0.113, 0.119, 0.12~

Part B

Use unique() to show all the variables present in your combined dataset.

unique(data$var)
##  [1] "GN_Al"        "GN_ANC"       "GN_Ca"        "GN_Cl"        "GN_DIC"      
##  [6] "GN_DOC"       "GN_H"         "GN_HCO3"      "GN_K"         "GN_Mg"       
## [11] "GN_Na"        "GN_NH4_N"     "GN_NO3_N"     "GN_pH"        "GN_Si"       
## [16] "GN_SO4_S"     "GN_spCond"    "GN_TN"        "GN_TP"        "IS_discharge"

Q3

Before we tackle our main question, let’s do some exploratory data analysis (or EDA). This is a common first step when trying to answer a practical or scientific question. Rather than rushing into applying an analysis, we want to get a gestalt understanding of our data.

For our EDA, we will be making some plots. Don’t worry too much about proper labels or titles, as the purpose of these is to inform you, the investigator.

The plots should still be easily to interpret.

Part A

Let’s say you don’t trust my summary of the experimental maniupulation effects. Plot a facet-wrapped time series of pH, total nitrogen, and sulfate at the both sites.

(Hint: The metadata for all Macrosheds data is available in the ‘data/macrosheds_vardata.csv’)

data %>%
 filter(var %in% c("GN_pH", "GN_SO4_S", "GN_TN")) %>%
  ggplot(., aes(x = datetime, y = val, color = site_code))+
    geom_point()+
    facet_wrap(~var, scales = 'free', ncol = 1)

Part B

You don’t know the researchers at BBWM (I assume), so let’s do a quick validation of their measurements.

Make a scatterplot of the ratio of bicarbonate to total dissolved inorganic carbon vs pH at the experimental watershed.

data %>%
  pivot_wider(id_cols = c(datetime, site_code),
              names_from = var,
              values_from = val) %>%
  select(datetime, site_code, ph = GN_pH, bicarb = GN_HCO3, dic = GN_DIC) %>%
  filter(site_code == 'WB') %>%
  ggplot(., aes(x = ph, y = bicarb/dic))+
  geom_point()
## Warning: Removed 2909 rows containing missing values (geom_point).

Part C

Does your graph from the previous part make you feel confident in this data? Why or why not? Cite figures or ideas from lectures in your answer.

Q4

Which watershed is weathering more?

To answer this question we need to approximate weathering rates at both watersheds. We’ll be calculating annual fluxes of silicon as a proxy for chemical weathering at both sites.

Part A

Why are we using silicon fluxes and not one of the other variables you listed in question 2, part B?

Part B

What is the frequency of the discharge data? What about the frequency of the silicon samples? Do they match? Prove it via code. Make sure your output is able to be understood by a person.

data %>%
  filter(var == 'IS_discharge') %>%
  mutate(diff = datetime-lag(datetime)) %>%
  select(var,diff) %>%
  tail()
## # A tibble: 6 x 2
##   var          diff  
##   <chr>        <drtn>
## 1 IS_discharge 1 days
## 2 IS_discharge 1 days
## 3 IS_discharge 1 days
## 4 IS_discharge 1 days
## 5 IS_discharge 1 days
## 6 IS_discharge 1 days
data %>%
  filter(var == 'GN_Si') %>%
  mutate(diff = datetime-lag(datetime)) %>%
  select(var,diff) %>%
  tail()
## # A tibble: 6 x 2
##   var   diff  
##   <chr> <drtn>
## 1 GN_Si 1 days
## 2 GN_Si 1 days
## 3 GN_Si 1 days
## 4 GN_Si 1 days
## 5 GN_Si 1 days
## 6 GN_Si 1 days
# Yes, they match at a daily interval.

Part C

Calculate the flux of silcon at the daily timescale both EB and WB. Filter your data to only WY1998-WY2002. Use glimpse() to show your result.

The area of WB is 10.06 ha and the area of EB is 10.76 ha.

Do not remove NAs! We will need those later!

(Hint: The metadata for all Macrosheds data is available in the ‘data/macrosheds_vardata.csv’)

flux <- data %>%
  pivot_wider(id_cols = c(datetime, site_code),
              names_from = var,
              values_from = val) %>%
  select(datetime, site_code, si = GN_Si, q = IS_discharge) %>%
  filter(datetime > as.POSIXct('1998-10-01'),
         datetime < as.POSIXct('2002-09-30')) %>%
  mutate(area = ifelse(site_code == 'EB', 10.76, 10.06),
         flux = ((si*q)/area)*(86400/1)*(1/1e6)) # calc flux in mg/sec*ha and convert sec -> day and mg -> kg

Part D

Plot your new, fine-scale flux timeseries with both line and point geometry on the same plot. Be sure to include appropriate labels, a title, and a descriptive caption. Color your timeseries to match the map in the background section.

ggplot(flux, aes(x = datetime, y = flux, color = site_code))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c('blue', 'red'))+
  labs(y = 'Si Flux (kg/ha/day)', x = 'Date', color = 'Site', title = 'Daily Silicon Flux at BBWM',
       caption = 'Daily fluxes of silicon, a rock derived solute, at BBWM. \n WB has been treated with acid, EB is a reference watershed.')
## Warning: Removed 1251 rows containing missing values (geom_point).
## Warning: Removed 232 row(s) containing missing values (geom_path).

Part E

Is our data completely continuous? What will happen if we try to calculate flux with this dataset? Is there a pattern to the breaks? Why might the good folks at BBWM let those breaks happen?

(Hint: Look at timeseries of the Q and Si separately. Think the flux formula and how you would spend your research dollars.)

Part F

Let’s fill in our to make it continuous. You should be careful about doing this in your own work, but for reasons you (should) have discussed in part E we can safely do it here. There are many ways to do (with varying degrees of merit). We will be using a simple one: linear interpolation. Essentially, we are telling R to draw a straight line between observations on either side of the data gap.

To do this, we will use the na.approx() function. The ‘rule’ for linear interpolation is 2.

flux_full <- flux %>%
  arrange(datetime) %>%
  mutate(flux_int = na.approx(flux, rule = 2))

Part G

Plot your new, interpolated timeseries just like you did for part D.

p <- ggplot(flux_full, aes(x = datetime, y = flux_int, color = site_code))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c('blue', 'red'))+
  labs(y = 'Si Flux (kg/ha/day)', x = 'Date', color = 'Site', title = 'Daily Silicon Flux at BBWM',
       caption = 'Daily fluxes of silicon, a rock derived solute, at BBWM. \n WB has been treated with acid, EB is a reference watershed.')

ggplotly(p)

Part H

Zoom in on the previous gaps of your graph. Using such a naive method is… not great. But it’s fine for our purposes and it shouldn’t throw off our annual estimates by much.

We finally have what we need to calculate annual silicon flux for both watersheds. Do so and present your data as you find appropriate in units of kg per ha per year. If you color your data, make sure it sticks to the schema we started earlier.

Your figure or table should clearly answer our question: which watershed hosts the higher weathering rates?

Hints:

  • If you are struggling check out the work with data primer on Rstudio cloud.

  • The function water_year() slickly pulls the WY from a date.

  • Use a caption to explictly state which watershed has higher weathering rates.

flux_full %>%
  mutate(wy = water_year(datetime)) %>%
  group_by(site_code, wy) %>%
  summarize(flux = sum(flux_int)) %>%
  ggplot(., aes(x = as.character(wy), y = flux, color = site_code))+
    geom_point()+
  scale_color_manual(values = c('blue', 'red'))+
  labs(y = 'Si Flux (kg/ha/year)', x = 'Date', color = 'Site', title = 'Annual Silicon Flux at BBWM',
       caption = 'Annual fluxes of silicon, a rock derived solute, at BBWM. \n WB has been treated with acid, EB is a reference watershed. \n Note the higher fluxes at the treatment watershed.')
## `summarise()` regrouping output by 'site_code' (override with `.groups` argument)

Part I

Consider your two time series (the one with gaps and the interpolated one) and your annual estimates. Do you think your estimates are biased high, low, or are perfect? Why do you think that?

Does this affect your conclusion from the previous part? Why or why not?

Bonus

Go take a look at the Gaillardet et al 1999 paper in the ‘pubs’ folder of this project. This paper is super dense, so don’t try to understand every bit of it. A good method for gleaming understanding from scary papers is to read the abstract, conclusions, and figures. Before we go any further, please give the paper a hearty skim.

Part A

One way Gaillardet and his co-authors understand solute fluxes in large rivers of the world is through mixing ratios. The basic idea is that differing parent materials (rocks that underlay soils) have different elemental compositions that in turn yield identifiable solute flux signatures.

Recreate figure 2 from the paper for BBWM, with each year of record at EB/WB as a point. Color your data by site and include a 1:1 reference line.

Part B

What do your mixing ratio plots indicate is the underlying parent material of BBWM?

Part C

What may be skewing our diagram and how could we correct it with addtional data?

(Hint: Think about where BBWM is and how mixing ratios are calculated. Are there any potential non-weathering sources around?)